** To activate the ICCs macro, you need to specify where the macro is located on
   your computer.  This example program specifies the location as a network
   drive (S), and the subdirectory path.
**;

*%include '<drive letter>:\<subdirectory path>\iccs.sas';

** To call the macro, once it has been activated, you need to identify a data
   set for analysis.  If the data is permanent, you will need to use a libname
   statement to point to the directory where the file is stored.
**;

*libname <name of library> '<path to data set>';

** For a temporary data set, already created during the current session, WORK
   is the appropriate library reference (this is the default library reference).
**;

** Example using data from Shrout and Fleiss **;

%include 's:\paul weiss\haber\iccs.sas' ;

data sandf;
   input id judge1 judge2 judge3 judge4;
   cards;
   1 9 2 5 8 
   2 6 1 3 2 
   3 8 4 6 8 
   4 7 1 2 6 
   5 10 5 6 9 
   6 6 2 4 7
   ;
run;

%ICCs (lib=work, dsn=sandf, rep=rep, var=judge, nummeas=4, numsub=6, numrep=1, 
       varlist=judge1 judge2 judge3 judge4, dataname=Shrout and Fleiss Data,
       seed=12345, alpha=0.05) ;

*************************************************************************************;
*************************************************************************************;
*************************************************************************************;
***                                                                               ***;
*** ICCs Macro (SAS Version 9.1) -- Most recent update, March 21, 2005            ***;
*** Written by -- Paul Weiss and George Cotsonis                                  ***;
*** Under the Direction of -- Michael Haber                                       ***;
*** Supported by NIH award R01 MH070028                                           ***;
***                                                                               ***;
*** This macro will generate measures of agreement and consistency between raters ***;
*** using moments and components of variance. The macro requires that the data be ***;
*** in columnar form, with each individual observation representing a case, and   ***;
*** each column representing the measurements taken by an individual rater.       ***;
*** Each case needs a unique identifier, called ID, which must appear in the      ***;
*** data set as a variable.                                                       ***;
***                                                                               ***;
***                                                                               ***;
*** Example:                                                                      ***;
***                                                                               ***;
*** ID Measure1 Measure2 Measure3 Measure4                                        ***;
***  1    8         8        9        9                                           ***;
***  2    9         9        8        8                                           ***;
***  3    8         9        9        8                                           ***;
***  4    9         8        8        9                                           ***;
***                                                                               ***;
*** The macro requires the user to specify the following parameters:              ***;
***                                                                               ***;
*** lib     : the library where the data are stored (WORK is default)             ***;
*** dsn     : the data set for analysis                                           ***;
*** rep     : replicate variable name [Not Working]                               ***;
*** var     : a one-word description of the measurement                           ***;
*** nummeas : the number of measurers                                             ***;
*** numsub  : the number of individuals being measured                            ***;
*** numrep  : the number of replicates on each individual (default is 1)          ***;
*** varlist : an explicit listing of variables for the analysis                   ***;
*** dataname: explicit description of data source                                 ***;
*** seed    : random number seed for method 2 (random selection of measure)       ***;
*** alpha   : significance level for confidence intervals                         ***;
***                                                                               ***;
*************************************************************************************;
*************************************************************************************;
*************************************************************************************;

%macro ICCs (lib=, dsn=, rep=, var=, numrater=, nummeas=, numsub=, numrep=, varlist=, dataname=, 
             seed=0, alpha=0.05 ) ;

footnote;

*** Page 1: First 25 observations in the data ***;

Title1 ;
title2 ;
title3 "&dataname -- &var";

%if &numsub gt 25 %then %do;
Title4 "First 25 Observations in the Data";
%end; %else %do;
Title4 "First &numsub Observations in the Data";
%end;

proc print data=&lib..&dsn (obs=25);
run;

%if &numrep eq 1 %then %do; *** Begin NON-REPLICATE ICC calculations ***;

*** Page 2: Means And Standard Deviations for each observer ***;

Title1 ;
title2 ;
title3 "&dataname -- &var";
Title4 "Simple Descriptive Statistics For Each Observer";

proc means mean std data=&lib..&dsn noprint;
   var &varlist;
   output out=means mean=
   %do bar=1 %to &nummeas;
   mean&var&bar
   %end;
   std=
   %do sd=1 %to &nummeas;
   sd&var&sd
   %end;
   ;
run;

%do ds=1 %to &nummeas;

data rater&ds;
   set means;
   mean=mean&var&ds;
   sd  =sd&var&ds;
   keep mean sd;
run;

%end;

data raters;
   set
   %do ds=1 %to &nummeas;
   rater&ds
   %end;
   ;
   rater+1;
   label rater="Observer"
         mean="Mean"
         sd="Standard Deviation"
   ;
run;

proc print noobs label;
   var rater mean sd;
run;


*** Start with Pearson Correlations ***;

Title1 ;
title2 ;
title3 "&dataname -- &var";
Title4 "Pearson Correlation Coefficients and Covariances";

*** Preparing data for analysis ***;

ods select none;
proc corr data=&lib..&dsn cov;
   var &varlist;
   ods output cov=covmat;
run;

%do d=1 %to &nummeas;

data a&d;
   set covmat;
   if _n_=&d;
   %do j=1 %to &nummeas;
   &var&d&j=&var&j;
   %end;
run;

%end;

data oneline;
   merge
   %do d=1 %to &nummeas;
   a&d
   %end;
   ;;
   keep
   %do i=1 %to &nummeas;
       %do j=1 %to &nummeas;
	       &var&i&j
       %end;
   %end;
   ;; 
run;

data pear;
   set oneline;

   %do i=1 %to &nummeas;
       %do j=1 %to &nummeas;

	   r&i&j = &var&i&j / sqrt(&var&i&i * &var&j&j);

	   %end;
   %end;

run;

%let a11 = 11;

data corrs;
   set pear;
   array r(&nummeas,&nummeas) r11 -- r&nummeas&nummeas;
   array x(&nummeas,&nummeas) &var&a11 -- &var&nummeas&nummeas;
   do i=1 to &nummeas;
      do j=1 to &nummeas;
	     pearson=r(i,j);
		 cov=x(i,j);
		 var1=x(i,i);
		 var2=x(j,j);
		 if i lt j then output;
	  end;
   end;
   keep i j pearson cov var1 var2;
   label pearson="Pearson Correlation"
         cov="Covariance"
         var1="Variance 1"
		 var2="Variance 2"
		 i="Measure 1"
		 j="Measure 2"
		 ;
run;

ods select all;
proc print label noobs;
   var i j pearson cov ;
run;

*** Calculating the OCCC ***;

Title4 "Overall Concordance Correlation Coefficients";

* Need Covariances from Pearson *;

proc means sum noprint data=corrs;
   var cov;
   output out=covs sum=sumcov;
run;

proc sort data=corrs;
   by i j;
run;

data v1;
   set corrs;
   by i j;
   if first.i; *** Keeps first j-1 variances ***;
   measure=i;
   var=var1;
   keep measure var;
run;
   
data v2;
   set corrs;
   by i j;
   if last.j; *** Keeps last j-1 variances ***;
   measure=j;
   var=var2;
   keep measure var;
run;

data v;
   set v1 v2;
run;

proc sort;
   by measure;
run;

data vx;
   set v;
   by measure;
   if first.measure;
run;

proc means sum noprint data=vx;
   var var;
   output out=vars sum=sumvar;
run;

proc means mean data=&lib..&dsn noprint;
   var &varlist;
   output out=means mean=;
run;

data meandiffs;
   set means;
   array means (&nummeas) &varlist;
   xend = &nummeas - 1;
   meansq=0; * to start *;
   do i = 1 to xend;
      do j = i+1 to &nummeas;
	     meansq + (means(i) - means(j))**2;
	  end;
   end;
   drop i j xend;
run;

proc means mean data=corrs noprint;
   var pearson;
   output out=meancorr mean=corrbar;
run;

data occc;
   merge covs vars meandiffs meancorr;
   nab = 2*sumcov / (((&nummeas - 1)*sumvar) + meansq);
   nad = 2*sumcov / (((&nummeas - 1)*sumvar));
   nl  = corrbar;
   label nab = "Coefficient of Absolute Agreement"
         nad = "Coefficient of Additive Agreement"
		 nl  = "Coefficient of Linear Agreement"
         sumcov = "Sum of Covariance"
		 sumvar = "Sum of Variance"
		 meansq  = "Sum of Squared Difference in Means"
		 ;
proc print label noobs;
   var sumcov sumvar meansq nab nad nl;
run;

*** Calculating ICC-C and ICC-A ***;

Title4 "Intraclass Correlation Coefficients for Agreement and Consistency";


data icc;

   set &lib..&dsn;

   %do i=1 %to &nummeas;

       id=id; &var=&var&i; rater=&i; output;

   %end;
   keep id &var rater;

run;

ods select none;
proc glm data=icc;
   class id rater ;
   model &var = id rater ;
   random id ;
   ods output OverallANOVA=over ModelANOVA=SS ;
run;

data one;
   set over;
   if Source="Error";
   edf=df; ems=ms;
   keep edf ems;
run;

data two;
   set ss;
   if Hypothesistype=3;
   if source='id' then do;
      bdf=df; bms=ms; jdf=0; jms=0; fobs=fvalue;
   end; else
   if source='rater' then do;
      bdf=0; bms=0; jdf=df; jms=ms; fobs=0;
   end;

proc means sum noprint data=two;
   var bdf bms jdf jms fobs;
   output out=twoa sum= ;
run;

data iccfull;
   merge one twoa;
   drop _type_ _freq_;

   n=&numrep;    * This will be assigned in a macro variable later as number of replicates *;
   a=&nummeas;   * This will be assigned in a macro variable later as number of raters     *;
   b=&numsub;    * This will be assigned in a macro variable later as number of subjects   *;

   sse = ems;

   sss = (bms - ems)/(n * a);

   ssr = (jms - ems)/(n * b);

   ICCa = sss / (sss + ssr + sse);
   ICCc = sss / (sss + sse);

   conf_lim=1-(&alpha/2);

   *** ICC-C Confidence Limits ***;
   cdfnum=b-1;
   cdfden=(b-1)*(a-1);

   cftabledl=finv(conf_lim, cdfnum, cdfden);
   cftabledu=finv(conf_lim, cdfden, cdfnum);

   fl = fobs/cftabledl;
   fu = fobs*cftabledu;

   ccll= (fl-1)/(fl + a - 1);
   cclu= (fu-1)/(fu + a - 1);

   *** Confidence Limits for ICC-A ***;

   aa=(a*ICCa)/(b*(1-ICCa));
   ab=1+((a*ICCa*(b-1))/(b*(1-ICCa)));

   asatnum=(aa*jms+ab*ems)**2;
   asatden1=((aa*jms)**2)/(a-1);
   asatden2=((ab*ems)**2)/((a-1)*(b-1));

   anu = asatnum/(asatden1 + asatden2);

   aftabledl=finv(conf_lim, a-1, anu);
   aftabledu=finv(conf_lim, anu, a-1);

   aclltop=b*(bms-aftabledl*ems);
   acllbot=aftabledl*(b*jms+(ab - a -b)*ems)+b*bms;

   acll=aclltop/acllbot;

   aclutop=b*(aftabledu*bms-ems);
   aclubot=b*jms+(ab - a -b)*ems + b*aftabledu*bms;

   aclu=aclutop/aclubot;

   rater1="Overall";
   rater2=" ";

   label sse="Error Component"
         ssr="Rater Component"
		 sss="Subject Component"
		 ICCa="Intraclass Correlation (Agreement)"
		 ICCc="Intraclass Correlation (Consistency)"
		 ccll="Lower Confidence Limit for ICC-C (alpha=&alpha)"
		 cclu="Upper Confidence Limit for ICC-C (alpha=&alpha)"
		 acll="Lower Confidence Limit for ICC-A (alpha=&alpha)"
		 aclu="Upper Confidence Limit for ICC-A (alpha=&alpha)"
   ;

run;


%if &nummeas ge 3 %then %do; 

%do x=1 %to %eval(&nummeas - 1);
%do y=%eval(&x+1) %to &nummeas;

data posthoc; *** Create a data set with only two measures, repeat for all pairwise measures ***;
   set &lib..&dsn;
   var1=&var&x; * Pick first measure *;
   var2=&var&y; * Pick second measure *;
   keep id var1 var2; * Keep data set with chosen measures as measure 1 and measure 2 *;
run;


data icc;

   set posthoc;

   %do i=1 %to 2;

       id=id; var=var&i; rater=&i; output;

   %end;
   keep id var rater;

run;

ods select none;
proc glm data=icc;
   class id rater ;
   model var = id rater ;
   random id ;
   ods output OverallANOVA=over ModelANOVA=SS ;
run;

data one;
   set over;
   if Source="Error";
   edf=df; ems=ms;
   keep edf ems;
run;

data two;
   set ss;
   if Hypothesistype=3;
   if source='id' then do;
      bdf=df; bms=ms; jdf=0; jms=0; fobs=fvalue;
   end; else
   if source='rater' then do;
      bdf=0; bms=0; jdf=df; jms=ms; fobs=0;
   end;

proc means sum noprint data=two;
   var bdf bms jdf jms fobs;
   output out=twoa sum= ;
run;

data iccpost&x&y;
   merge one twoa;
   drop _type_ _freq_;

   n=&numrep;    * This will be assigned in a macro variable later as number of replicates *;
   a=2       ;   * Use two (2) for pairwise post-hoc analysis in ICC calculations          *;
   b=&numsub;    * This will be assigned in a macro variable later as number of subjects   *;

   sse = ems;

   sss = (bms - ems)/(n * a);

   ssr = (jms - ems)/(n * b);

   ICCa = sss / (sss + ssr + sse);
   ICCc = sss / (sss + sse);

   conf_lim=1-(&alpha/2);

   *** ICC-C Confidence Limits ***;
   cdfnum=b-1;
   cdfden=(b-1)*(a-1);

   cftabledl=finv(conf_lim, cdfnum, cdfden);
   cftabledu=finv(conf_lim, cdfden, cdfnum);

   fl = fobs/cftabledl;
   fu = fobs*cftabledu;

   ccll= (fl-1)/(fl + a - 1);
   cclu= (fu-1)/(fu + a - 1);

   *** Confidence Limits for ICC-A ***;

   aa=(a*ICCa)/(b*(1-ICCa));
   ab=1+((a*ICCa*(b-1))/(b*(1-ICCa)));

   asatnum=(aa*jms+ab*ems)**2;
   asatden1=((aa*jms)**2)/(a-1);
   asatden2=((ab*ems)**2)/((a-1)*(b-1));

   anu = asatnum/(asatden1 + asatden2);

   aftabledl=finv(conf_lim, a-1, anu);
   aftabledu=finv(conf_lim, anu, a-1);

   aclltop=b*(bms-aftabledl*ems);
   acllbot=aftabledl*(b*jms+(ab - a -b)*ems)+b*bms;

   acll=aclltop/acllbot;

   aclutop=b*(aftabledu*bms-ems);
   aclubot=b*jms+(ab - a -b)*ems + b*aftabledu*bms;

   aclu=aclutop/aclubot;

   rater1="&x";
   rater2="&y";

   label sse="Error Component"
         ssr="Rater Component"
		 sss="Subject Component"
		 ICCa="Intraclass Correlation (Agreement)"
		 ICCc="Intraclass Correlation (Consistency)"
		 ccll="Lower Confidence Limit for ICC-C (alpha=&alpha)"
		 cclu="Upper Confidence Limit for ICC-C (alpha=&alpha)"
		 acll="Lower Confidence Limit for ICC-A (alpha=&alpha)"
		 aclu="Upper Confidence Limit for ICC-A (alpha=&alpha)"
		 rater1="Rater 1"
		 rater2="Rater 2"
   ;

run;
%end;
%end;
%end;

data icc;
   set iccfull 
   %if &nummeas ge 3 %then %do;
      %do x=1 %to %eval(&nummeas - 1);
          %do y=%eval(&x+1) %to &nummeas;
   iccpost&x&y
          %end;
      %end;
   %end;
   ;;
run;

ods select all;
title5;

proc print label noobs data=icc;
   var %if &nummeas ge 3 %then %do;
       rater1 rater2 sss ssr sse ;
       %end; %else %do;
       sss ssr sse ;
	   %end;
	   ;;
run;

proc print label noobs data=icc;
   var %if &nummeas ge 3 %then %do;
       rater1 rater2 icca acll aclu ;
       %end; %else %do;
       icca acll aclu ;
	   %end;
	   ;;
run;

proc print label noobs data=icc;
   var %if &nummeas ge 3 %then %do;
       rater1 rater2 iccc ccll cclu;
       %end; %else %do;
       iccc ccll cclu;
	   %end;
	   ;;
run;

%end;  * End NON-REPLICATE Coefficient calculation *;

%if &numrep ge 2 %then %do; *** Only do this if there are replicates in the data ***;

**** Begin block for calculating CIV ****;

Title4 "Coefficient of Interobserver Variability from Sums of Squares";

data CIV;
    set &lib..&dsn;

    **** set up arrays, y(j,k), yj=y(j,.), u(j);
    array y(&numrater,&numrep) &varlist;
    array yj(&numrater) yj1-yj&numrater;
    array u(&numrater) u1-u&numrater;
 
    ******** initialize ydd ( y..);

    ydd=0;

    ********* calculate v;
    do j=1 to &numrater;
	   yj(j)=0;  ******* this is yj.;
	   do k=1 to &numrep;
	      ydd=ydd+y(j,k);   ***** get ydd.;
		  yj(j)=yj(j)+y(j,k);  **** get yj.;
       end;
	end;

	*** Calculate the means for ydd and yj ***;

	yddbar = ydd/%eval(&numrater*&numrep);

	array yjbar (&numrater) yjbar1-yjbar&numrater ;

	do j=1 to &numrater;
	   yjbar(j)=yj(j)/&numrep ;
	end;

	v=0;

    do j=1 to &numrater;
	   v=v+(yjbar(j)-yddbar)**2;
    end;

    v=v/(&numrater-1);

    *****  calculate u(j) and u. ;
	**** u. is used in calculation of t2 and w2;
    ud=0;

    do j=1 to &numrater;
	   u(j)=0;
	   do k=1 to &numrep;
	      u(j)=u(j)+(y(j,k)-yjbar(j))**2;
	   end;
	   u(j)=u(j)/(&numrep-1);
	   ud=ud+u(j);
	end;

	ud=ud/(&numrater);

    keep v ud u1 - u&numrater;

run;

proc means noprint;
   var v ud u1 - u&numrater;
   output out=CIVa mean=vd udd u1 - u&numrater;  **** all we need are means;
run;

data CIVb;
   set CIVa;
   t2=vd-udd/&numrep;
   w2=udd;
   civ=t2/(t2+w2);
   psi=1-civ;
   keep vd udd t2 w2 civ psi u1 - u&numrater;
   label vd='V.'
	     udd='U..'
		 t2='Tau Squared'
		 w2='Omega Squared'
		 civ='CIV'
		 psi='PSI'
   ;
run;

proc print label noobs;
   var vd udd t2 w2 civ psi;
run;


**** Moments and CCCs, Method 1 ****;

Title4 "Moments and CCCs, Method 1";
title5 "Based on the First Observation of each Observer";

data mccc1;

   set &lib..&dsn;

   array y (&numrater, &numrep) &varlist ;
   array frst (&numrater) f1-f&numrater ;

   do j=1 to &numrater;
      frst(j)=y(j,1);
   end;

   keep f1 - f&numrater;

run;

proc means noprint data=mccc1;
   var f1 - f&numrater;
   output out=mccc1a mean=m1 - m&numrater
                     var =v1 - v&numrater
   ;
run;

data mccc1b;
   merge mccc1a CIVb;

   array sigma (&numrater) v1 - v&numrater ;
   array m (&numrater) m1 - m&numrater ;
   array w (&numrater) w1 - w&numrater ;
   array d (&numrater) d1 - d&numrater ;
   array u (&numrater) u1 - u&numrater ;

   do j=1 to &numrater;

      w(j) = u(j);
	  d(j) = sigma(j) - w(j);

   end;

run;

data mccc1c;
   set mccc1b;

   array sigma (&numrater) v1 - v&numrater ;
   array m (&numrater) m1 - m&numrater ;
   array w (&numrater) w1 - w&numrater ;
   array d (&numrater) d1 - d&numrater ;
   array u (&numrater) u1 - u&numrater ;

   do j=1 to &numrater ;

      rater=j;
      muhat=m(j);
	  deltahatsquared=d(j);
	  whatsquared=w(j);
	  sigmahatsquared=sigma(j);
	  sigmahat=sqrt(sigma(j));

	  output;

   end;

   label rater="Rater"
         muhat="Mu"
		 deltahatsquared="Delta Squared"
		 whatsquared="Omega Squared"
		 sigmahatsquared="Sigma Squared"
		 sigmahat="Sigma"
   ;

run;

proc print noobs label;
   var rater muhat deltahatsquared whatsquared sigmahatsquared sigmahat;
run;
 
ods select none; *** Want to keep CORR matrix, but not output to screen ***;
proc corr data=mccc1 ;
   var f1 - f&numrater;
   ods output PearsonCorr=corrs;
run;

%do d=1 %to &numrater;

data r&d;
   set corrs;
   if _n_=&d;
   %do j=1 %to &numrater;
   r&d&j=f&j;
   %end;
run;

%end;

data r1line;
   merge
   %do d=1 %to &numrater;
   r&d
   %end;
   ;;
   keep
   %do i=1 %to &numrater;
       %do j=1 %to &numrater;
	       r&i&j
       %end;
   %end;
   ;;

run;

data mccc1d;
   merge mccc1b r1line;

   array sigma (&numrater) v1 - v&numrater ;
   array m (&numrater) m1 - m&numrater ;
   array w (&numrater) w1 - w&numrater ;
   array d (&numrater) d1 - d&numrater ;
   array u (&numrater) u1 - u&numrater ;

   array rho (&numrater,&numrater) r11 -- r&numrater&numrater ;

   do j=1 to %eval(&numrater-1) ;
      do jprime=j+1 to &numrater ;

	     rater1 = j;
		 rater2 = jprime;

		 rhohat=rho(j,jprime);

		 rhohatmu = rho(j,jprime) * sqrt( sigma(j)*sigma(jprime) / (d(j)*d(jprime)) ) ;

		 output;

	  end;
   end;


   label rater1="Rater 1"
         rater2="Rater 2"
		 rhohat="Pearson Correlation of First Measures"
		 rhohatmu="Rho-mu"
   ;
   
run;

data mccc1f;
   merge mccc1b r1line;

   array sigma (&numrater) v1 - v&numrater ;
   array m (&numrater) m1 - m&numrater ;
   array w (&numrater) w1 - w&numrater ;
   array d (&numrater) d1 - d&numrater ;
   array u (&numrater) u1 - u&numrater ;
   array rho (&numrater,&numrater) r11 -- r&numrater&numrater ;

   denompij=0;

   do j=1 to %eval(&numrater-1);
      do jprime = j+1 to &numrater;
	     denompij=denompij + (rho(j,jprime) * sqrt(sigma(j) * sigma(jprime)) / 
                             (&numrater * (&numrater-1)/2));
	  end;
   end;

   do j=1 to &numrater;
      rater=j;
	  mu_j=m(j);
	  delta_sq_j = d(j);
	  omega_sq_j = w(j);
	  sigma_sq_j = sigma(j);
      rhohatij= d(j)/sigma(j) ;
	  pij = sigma(j) / denompij;
	  checker=pij*(1-rhohatij);
	  output;
   end;

   label mu_j="Mu"
         delta_sq_j="Delta- Squared"
         omega_sq_j="Omega- Squared"
         sigma_sq_j="Sigma- Squared"
         rhohatij="Rho^I"
		 pij="Pi"
		 checker="Product for checking rho-hats"
		 ;

run;



**** Overall Measures of Agreement, Method 1 ****;
**** Calculated using first replicate        ****;

**** Use code from non-rep (above) and modify to handle first obs from data ****;

data cccr;
   set &lib..&dsn;
   array y (&numrater, &numrep) &varlist ;
   array frst (&numrater) f1-f&numrater ;

   do j=1 to &numrater;
      frst(j)=y(j,1);
   end;

   keep f1 - f&numrater;

run;

ods select none;
proc corr cov;
   var f1-f&numrater;
   ods output cov=covmat;
run;

%do d=1 %to &numrater;

data a&d;
   set covmat;
   if _n_=&d;
   %do j=1 %to &numrater;
   &var&d&j=f&j;
   %end;
run;

%end;

data oneline;
   merge
   %do d=1 %to &numrater;
   a&d
   %end;
   ;;
   keep
   %do i=1 %to &numrater;
       %do j=1 %to &numrater;
	       &var&i&j
       %end;
   %end;
   ;; 
run;

data pear;
   set oneline;

   %do i=1 %to &numrater;
       %do j=1 %to &numrater;

	   r&i&j = &var&i&j / sqrt(&var&i&i * &var&j&j);

	   %end;
   %end;

run;

%let a11 = 11;

data corrs;
   set pear;
   array r(&numrater,&numrater) r11 -- r&numrater&numrater;
   array x(&numrater,&numrater) &var&a11 -- &var&numrater&numrater;
   do i=1 to &numrater;
      do j=1 to &numrater;
	     pearson=r(i,j);
		 cov=x(i,j);
		 var1=x(i,i);
		 var2=x(j,j);
		 if i lt j then output;
	  end;
   end;
   keep i j pearson cov var1 var2;
   label pearson="Pearson Correlation"
         cov="Covariance"
         var1="Variance 1"
		 var2="Variance 2"
		 i="Measure 1"
		 j="Measure 2"
		 ;
run;

*** Calculating the OCCC ***;

* Need Covariances from Pearson *;

proc means sum noprint data=corrs;
   var cov;
   output out=covs sum=sumcov;
run;

proc sort data=corrs;
   by i j;
run;

data v1;
   set corrs;
   by i j;
   if first.i; *** Keeps first j-1 variances ***;
   measure=i;
   var=var1;
   keep measure var;
run;
   
data v2;
   set corrs;
   by i j;
   if last.j; *** Keeps last j-1 variances ***;
   measure=j;
   var=var2;
   keep measure var;
run;

data v;
   set v1 v2;
run;

proc sort;
   by measure;
run;

data vx;
   set v;
   by measure;
   if first.measure;
run;

proc means sum noprint data=vx;
   var var;
   output out=vars sum=sumvar;
run;

proc means mean data=cccr noprint;
   var f1-f&numrater;
   output out=means mean=;
run;

data meandiffs;
   set means;
   array means (&numrater) f1-f&numrater;
   xend = &numrater - 1;
   meansq=0; * to start *;
   do i = 1 to xend;
      do j = i+1 to &numrater;
	     meansq + (means(i) - means(j))**2;
	  end;
   end;
   drop i j xend;
run;

proc means mean data=corrs noprint;
   var pearson;
   output out=meancorr mean=corrbar;
run;

data occc;
   merge covs vars meandiffs meancorr;
   nab = 2*sumcov / (((&numrater - 1)*sumvar) +  meansq);
   label nab = "Coefficient of Absolute Agreement"
		 ;
run;


data mccc1e;
   merge mccc1b r1line occc;

   array sigma (&numrater) v1 - v&numrater ;
   array m (&numrater) m1 - m&numrater ;
   array w (&numrater) w1 - w&numrater ;
   array d (&numrater) d1 - d&numrater ;
   array u (&numrater) u1 - u&numrater ;

   array rho (&numrater,&numrater) r11 -- r&numrater&numrater ;
   num=0;

   rhohatmuc=0;
   muss=(
   %do k = 1 %to %eval(&numrater-1);
   m&k +
   %end;
   m&numrater)/&numrater;
   do j=1 to %eval(&numrater-1) ;
      do jprime=j+1 to &numrater ;
         num=num+2*(sqrt((sigma(j)*sigma(jprime))*rho(j,jprime)));
	  end;
   end;
   denom=((&numrater*(
		 %do i=1 %to %eval(&numrater-1);
         ((m&i - muss)**2) +
		 %end;
		 (m&numrater - muss)**2) +
		 (&numrater-1)*(
		 %do i=1 %to %eval(&numrater-1);
         d&i +
		 %end;
         d&numrater) ));
   rhohatmuc=num/denom;
   
   label rhohatmuc="Rho-c (mu)"
         nab="Rho-c (y)"
		 ;    

run;

*** Check sum ***;

data check;
   set mccc1f;
   cs+(pij * (1 - rhohatIJ)/&numrater);
   if _n_ = &numrater;
   keep cs;
run;
   
data check2;
   merge mccc1e check;
   left=1/nab;
   right=(1/rhohatmuc) +cs;
   keep left right;
run;

**** Now we have a data set with the Ordinary CCCs ****;

ods select all;

proc print noobs label data=mccc1d;
   var rater1 rater2 rhohat rhohatmu ;
run;

proc print noobs data=mccc1f label;
   var rater mu_j delta_sq_j omega_sq_j sigma_sq_j pij rhohatij checker;
run;

proc print noobs label data=mccc1d;
   var rater1 rater2 rhohat rhohatmu ;
run;

proc print data=mccc1e label noobs;
   var nab rhohatmuc;
run;

proc print data=check2 noobs;
   var left right;
   title6 "Checking Rho-c(y) and Rho-c(mu)";
run;


ods select none;


**** END METHOD 1 BLOCK ****;






**** Moments and CCCs, Method 2 ****;

Title4 "Moments and CCCs, Method 2";
title5 "Based on a Randomly Selected Observation from each Observer";

data mccc1;

   set &lib..&dsn;

   array y (&numrater, &numrep) &varlist ;
   array frst (&numrater) f1-f&numrater ;

   do j=1 to &numrater;
      rand=rantbl(&seed,
	  %do k=1 %to (&numrep-1);
	  1/&numrep,
	  %end;
      1/&numrep);
      frst(j)=y(j,rand);
   end;

   *keep f1 - f&numrater;

run;

proc means noprint;
   var f1 - f&numrater;
   output out=mccc1a mean=m1 - m&numrater
                     var =v1 - v&numrater
   ;
run;

data mccc1b;
   merge mccc1a CIVb;

   array sigma (&numrater) v1 - v&numrater ;
   array m (&numrater) m1 - m&numrater ;
   array w (&numrater) w1 - w&numrater ;
   array d (&numrater) d1 - d&numrater ;
   array u (&numrater) u1 - u&numrater ;

   do j=1 to &numrater;

      w(j) = u(j);
	  d(j) = sigma(j) - w(j);

   end;

run;

data mccc1c;
   set mccc1b;

   array sigma (&numrater) v1 - v&numrater ;
   array m (&numrater) m1 - m&numrater ;
   array w (&numrater) w1 - w&numrater ;
   array d (&numrater) d1 - d&numrater ;
   array u (&numrater) u1 - u&numrater ;

   do j=1 to &numrater ;

      rater=j;
      muhat=m(j);
	  deltahatsquared=d(j);
	  whatsquared=w(j);
	  sigmahatsquared=sigma(j);
	  sigmahat=sqrt(sigma(j));

	  output;

   end;

   label rater="Rater"
         muhat="Mu"
		 deltahatsquared="Delta Squared"
		 whatsquared="Omega Squared"
		 sigmahatsquared="Sigma Squared"
		 sigmahat="Sigma"
   ;

run;

ods select all;

proc print noobs label;
   var rater muhat deltahatsquared whatsquared sigmahatsquared sigmahat;

run;
 
ods select none; *** Want to keep CORR matrix, but not output to screen ***;
proc corr data=mccc1 ;
   var f1 - f&numrater;
   ods output PearsonCorr=corrs;
run;

%do d=1 %to &numrater;

data r&d;
   set corrs;
   if _n_=&d;
   %do j=1 %to &numrater;
   r&d&j=f&j;
   %end;
run;

%end;

data r1line;
   merge
   %do d=1 %to &numrater;
   r&d
   %end;
   ;;
   keep
   %do i=1 %to &numrater;
       %do j=1 %to &numrater;
	       r&i&j
       %end;
   %end;
   ;;

run;

data mccc1d;
   merge mccc1b r1line;

   array sigma (&numrater) v1 - v&numrater ;
   array m (&numrater) m1 - m&numrater ;
   array w (&numrater) w1 - w&numrater ;
   array d (&numrater) d1 - d&numrater ;
   array u (&numrater) u1 - u&numrater ;

   array rho (&numrater,&numrater) r11 -- r&numrater&numrater ;

   do j=1 to %eval(&numrater-1) ;
      do jprime=j+1 to &numrater ;

	     rater1 = j;
		 rater2 = jprime;

		 rhohat=rho(j,jprime);

		 rhohatmu = rho(j,jprime)*sqrt(sigma(j)*sigma(jprime)/(d(j)*d(jprime))) ;

		 output;

	  end;
   end;

   label rater1="Rater 1"
         rater2="Rater 2"
		 rhohat="Pearson Correlation of First Measures"
		 rhohatmu="Rho-mu"
   ;

run;

data mccc1f;
   merge mccc1b r1line;

   array sigma (&numrater) v1 - v&numrater ;
   array m (&numrater) m1 - m&numrater ;
   array w (&numrater) w1 - w&numrater ;
   array d (&numrater) d1 - d&numrater ;
   array u (&numrater) u1 - u&numrater ;
   array rho (&numrater,&numrater) r11 -- r&numrater&numrater ;

   denompij=0;

   do j=1 to %eval(&numrater-1);
      do jprime = j+1 to &numrater;
	     denompij=denompij + (rho(j,jprime) * sqrt(sigma(j) * sigma(jprime)) / 
                             (&numrater * (&numrater-1)/2));
	  end;
   end;

   do j=1 to &numrater;
      rater=j;
	  mu_j=m(j);
	  delta_sq_j = d(j);
	  omega_sq_j = w(j);
	  sigma_sq_j = sigma(j);
      rhohatij= d(j)/sigma(j) ;
	  pij = sigma(j) / denompij;
	  checker=pij*(1-rhohatij);
	  output;
   end;

   label mu_j="Mu"
         delta_sq_j="Delta- Squared"
         omega_sq_j="Omega- Squared"
         sigma_sq_j="Sigma- Squared"
         rhohatij="Rho^I"
		 pij="Pi"
		 checker="Product for checking rho-hats"
		 ;

run;


data cccr;
   set &lib..&dsn;
   array y (&numrater, &numrep) &varlist ;
   array frst (&numrater) f1-f&numrater ;

   do j=1 to &numrater;
      frst(j)=y(j,1);
   end;

   keep f1 - f&numrater;

run;

ods select none;
proc corr cov;
   var f1-f&numrater;
   ods output cov=covmat;
run;

%do d=1 %to &numrater;

data a&d;
   set covmat;
   if _n_=&d;
   %do j=1 %to &numrater;
   &var&d&j=f&j;
   %end;
run;

%end;

data oneline;
   merge
   %do d=1 %to &numrater;
   a&d
   %end;
   ;;
   keep
   %do i=1 %to &numrater;
       %do j=1 %to &numrater;
	       &var&i&j
       %end;
   %end;
   ;; 
run;

data pear;
   set oneline;

   %do i=1 %to &numrater;
       %do j=1 %to &numrater;

	   r&i&j = &var&i&j / sqrt(&var&i&i * &var&j&j);

	   %end;
   %end;

run;

%let a11 = 11;

data corrs;
   set pear;
   array r(&numrater,&numrater) r11 -- r&numrater&numrater;
   array x(&numrater,&numrater) &var&a11 -- &var&numrater&numrater;
   do i=1 to &numrater;
      do j=1 to &numrater;
	     pearson=r(i,j);
		 cov=x(i,j);
		 var1=x(i,i);
		 var2=x(j,j);
		 if i lt j then output;
	  end;
   end;
   keep i j pearson cov var1 var2;
   label pearson="Pearson Correlation"
         cov="Covariance"
         var1="Variance 1"
		 var2="Variance 2"
		 i="Measure 1"
		 j="Measure 2"
		 ;
run;

*** Calculating the OCCC ***;

* Need Covariances from Pearson *;

proc means sum noprint data=corrs;
   var cov;
   output out=covs sum=sumcov;
run;

proc sort data=corrs;
   by i j;
run;

data v1;
   set corrs;
   by i j;
   if first.i; *** Keeps first j-1 variances ***;
   measure=i;
   var=var1;
   keep measure var;
run;
   
data v2;
   set corrs;
   by i j;
   if last.j; *** Keeps last j-1 variances ***;
   measure=j;
   var=var2;
   keep measure var;
run;

data v;
   set v1 v2;
run;

proc sort;
   by measure;
run;

data vx;
   set v;
   by measure;
   if first.measure;
run;

proc means sum noprint data=vx;
   var var;
   output out=vars sum=sumvar;
run;

proc means mean data=cccr noprint;
   var f1-f&numrater;
   output out=means mean=;
run;

data meandiffs;
   set means;
   array means (&numrater) f1-f&numrater;
   xend = &numrater - 1;
   meansq=0; * to start *;
   do i = 1 to xend;
      do j = i+1 to &numrater;
	     meansq + (means(i) - means(j))**2;
	  end;
   end;
   drop i j xend;
run;

proc means mean data=corrs noprint;
   var pearson;
   output out=meancorr mean=corrbar;
run;

data occc;
   merge covs vars meandiffs meancorr;
   nab = 2*sumcov / (((&numrater - 1)*sumvar) + meansq);
   label nab = "Coefficient of Absolute Agreement"
		 ;
run;


data mccc1e;
   merge mccc1b r1line occc;

   array sigma (&numrater) v1 - v&numrater ;
   array m (&numrater) m1 - m&numrater ;
   array w (&numrater) w1 - w&numrater ;
   array d (&numrater) d1 - d&numrater ;
   array u (&numrater) u1 - u&numrater ;

   array rho (&numrater,&numrater) r11 -- r&numrater&numrater ;
   num=0;

   rhohatmuc=0;
   muss=(
   %do k = 1 %to %eval(&numrater-1);
   m&k +
   %end;
   m&numrater)/&numrater;
   do j=1 to %eval(&numrater-1) ;
      do jprime=j+1 to &numrater ;
         num=num+2*(sqrt((sigma(j)*sigma(jprime))*rho(j,jprime)));
	  end;
   end;
   denom=((&numrater*(
		 %do i=1 %to %eval(&numrater-1);
         ((m&i - muss)**2) +
		 %end;
		 (m&numrater - muss)**2) +
		 (&numrater-1)*(
		 %do i=1 %to %eval(&numrater-1);
         d&i +
		 %end;
         d&numrater) ));
   rhohatmuc=num/denom;
   
   label rhohatmuc="Rho-c (mu)"
         nab="Rho-c (y)"
		 ;    

run;

data check;
   set mccc1f;
   cs+pij * (1 - rhohatIJ)/&numrater;
   if _n_ = &numrater;
run;
   
data check2;
   merge mccc1e check;
   left=1/nab;
   right=1/rhohatmuc+cs;
   keep left right;
run;

ods select all;

proc print noobs label data=mccc1d;
   var rater1 rater2 rhohat rhohatmu ;
run;

proc print noobs data=mccc1f label;
   var rater mu_j delta_sq_j omega_sq_j sigma_sq_j pij rhohatij checker;
run;

proc print noobs label data=mccc1d;
   var rater1 rater2 rhohat rhohatmu ;
run;

proc print data=mccc1e label noobs;
   var nab rhohatmuc;
run;

proc print data=check2 noobs;
   var left right;
   title6 "Checking Rho-c(y) and Rho-c(mu)";
run;

ods select none;


**** END METHOD 2 BLOCK ****;





**** Overall Measures of Agreement from ANOVA model ****;

title4 "Measures of Agreement based on 2-way Mixed Model with Interaction";

data anovaprep;
   set &lib..&dsn;

   %do j=1 %to &numrater; *** repeat over raters ***;

       %do k=1 %to &numrep; *** repeat over replicates ***;

	       &var = &var&j&k; rater=&j; rep=&k; output; *** create a dataset with each measure on its
		                                              *** own line. ;

	   %end; *** replicate loop ***;

   %end; *** rater loop ***;

   keep rater subject y rep;

run;
 
      ************** anova method starts here;
ods select none;

proc glm data=anovaprep;
   class subject rater;
   model y=rater subject subject*rater  ;
   random subject / test;
   ods output modelanova=modelanova;
run;

data anovaout;
   set modelanova;
   if Hypothesistype=3;
run;

data sr;
   set anovaout;
   if source='rater' ;
   msr=ms;
   keep msr;
run;

data ss;
    set anovaout;
	if source='subject' ;
	mss=ms;
	keep mss;
run;

data ssr;
    set anovaout;
	if source='subject*rater' ;
	mssr=ms;
	mse=ms/fvalue;
	keep mssr mse;
run;

data combo;
   merge sr ss ssr; 
   ss=(mss-mse)/(&numrep*&numrater);
   sr=(msr-mssr)/(&numrep*&numsub);
   ssr=(mssr-mse)/&numrep;
   se=mse;
   civ=(sr+ssr)/(sr+ssr+se);
   rhoy=(ss-ssr/(&numrater-1))/(ss+sr+ssr+se);
   rhou=(ss-ssr/(&numrater-1))/(ss+sr+ssr);
   keep ss sr ssr se civ rhoy rhou;
   label civ='CIV'
         rhoy='Rho-c (y) From Variance Components'
         rhou='Rho-c (mu) From Variance Components'
		 ss="Subject Component"
		 sr="Rater Component"
		 ssr="Subject- Rater Interaction Component"
		 se="Error Component"
   ;
run;

ods select all;
proc print noobs label;
   var ss sr ssr se civ rhoy rhou;
run;

%end;  *** End block for replicated measures ***;

%mend ICCs;